Model Energy Landscapes of Low- Temperature Fluids: Dipolar Hard 

Spheres 
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An analytical model of non-Gaussian energy landscape of low-temperature fluids is developed 
based on the thermodynamics of the fluid of dipolar hard spheres. The entire excitation profile of 
the liquid, from the high temperatures to the point of ideal-glass transition, has been obtained from 
the Monte Carlo simulations. The fluid of dipolar hard spheres loses stability when reaching the 
point of ideal-glass transition transforming via a first-order transition into a columnar liquid phase 
of dipolar chains locally arranged in a body-centered tetragonal order. 
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Thermodynamic and dynamic properties of super- 
cooled liquids are often related to the excess of config- 
urations they possess relative to the crystalline phase [l| • 
The number of states of the liquid state is reflected by 
the configurational entropy S C (T) conventionally defined 
0, H, 13 as the logarithm of the density of states Q(E) 
taken at the energy equal to the internal energy E(T) at 
temperature T: 



S c (T)=]n [n(£7(T))]. 



(1) 



The energy integral of the density of states and the Boltz- 
mann factor provide the canonical configuration integral 



',{(3) = J n(E)e- pE dE, 



(2) 



where (3 is the inverse temperature. 

Stillinger [f| Q gave an alternative definition of the 
configurational entropy in terms of the density of states 
Qcfi(E) of inherent structures, i.e. minima of the total en- 
ergy of the system as a function of its all translational an 
rotational degrees of freedom (energy landscape). The 
canonical configuration integral is now given by integrat- 
ing over the basin depths <f> and the free energy of inter- 
basin vibrations F v ((b), which is a weak, approximately 
linear, function of <j> @, HI 



(3) 



Because of the central role of the configurational en- 
tropy in early theories of viscous liquids by Adam, Gibbs, 
and DiMarzio [l(| , and more recent random first-order 
transition models by Wolynes and co-workers [TTf| , much 
effort has been invested in calculating the Stillinger con- 
figurational entropy, S%(T) = In \fl<h(<l)( T))] , from com- 
puter simulations of model fluids [3, |g, 02 • A general 
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connection between S C (T) and Sf(T) is, however, still 
unclear. 

Analytical modeling of the landscape of supercooled 
liquids is to a large extent based on the ideas advanced 
for spin glasses with quenched disorder (l3l . [T^ . The 
most prominent role in application to structural glasses 
is played by Derrida's random-energy, or Gaussian land- 
scape, model (REM) [Tij . Computer simulations, mostly 
limited to temperatures above the mode-coupling criti- 
cal temperature, basically support the REM, in particu- 
lar the prediction of the 1/T falloff of the average basin 
energy </>(T) from its high-temperature plateau @, 0. 
However, general theoretical arguments p, [l6| and re- 
cent simulations [l2| suggest that the low-energy portion 
of the Gaussian landscape might be inaccurate. In par- 
ticular, combinatorial arguments suggest that the deriva- 
tive of the enumeration function, a t p((j>) = N^ 1 ln^^ (</>)], 
approaches infinity at the point of ideal glass transi- 
tion when the system runs out of configurations and 
(t<I>(4>ig) — (N is the number of liquid particles). This 
infinite derivative eliminates the ideal glass transition at 
a positive temperature Q. A phenomenological descrip- 
tion, patching together logarithmic and Gaussian enu- 
meration functions, was suggested to provide a correct 
low-energy portion of (70 (0) [17| . 

Generally, unlike the case of spin glasses, there has 
been a lack of simple, solvable models of structural 
glasses. In this Letter we propose a new model of non- 
Gaussian landscape of low-temperature fluids with no 
quenched disorder. Our derivation is based on the estab- 
lished thermodynamics of the model fluid of dipolar hard 
spheres (DHS) [l8]. This monoatomic glass former (in 
contrast to binary mixtures used in many recent simula- 
tions 0,111) is l° n g known to resist phase transformations. 
This fluid lacks the liquid-vapor transition decomposing 
instead into low-density dipolar chains Upon 
cooling, dipolar hard spheres transform into a ferroelec- 
tric fluid [21| , but the ferroelectric phase is stabilized by 
tin-coil boundary conditions used in simulations and can 
be prevented by using a low dielectric constant for the 
reaction-field or Ewald corrections for the cutoff of dipo- 
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lar interactions [22J . This strategy has been employed in 
this study. We have carried out the Monte Carlo (MC) 
simulations of the DHS fluid within the standard NVT 
Metropolis protocol and the reaction-field cutoff correc- 
tions. The reaction- field dielectric constant of crf = 10, 
below the lowest value ~ 18 permitting the ferroelectric 
phase has allowed us to eliminate the transition to 
liquid ferroelectric. The initial configuration was set up 
as a face-centered cubic lattice of 108 dipolar spheres. 
The cubic simulation box also helps to suppress crystal- 
lization of dipoles which do not favor highly symmetric 
lattice structures [23|, [24[. The data were collected for 
(2 - 10) x 10 7 MC cycles. 

Apart from avoiding crystallization, a significant ad- 
vantage of the DHS fluid is the existence of a simple 
analytical form for the free energy (3F(p) = — ln[Z (/?)]. 
The Pade truncation of perturbation series suggested by 
Stell and co-workers [25| turned out to be very success- 
ful when tested against simulations [18j . The free energy 
of dipolar hard spheres depends on two parameters, the 
reduced density p* — pa 3 and the reduced temperature 
T = k^Ta 3 /m 2 . Here, p is the number density, a is the 
hard-sphere diameter, and m is the dipole moment. All 
calculations and simulations here have been done at con- 
stant volume with p* = 0.8 thus reducing the number of 
variables to one. 

Stell's Pade solution for F(/3) is 



(3F(J3) = Nf(rj) 



A/3 2 
1 + 6/3' 



(4) 



Here f(rj) is the reduced free energy of the fluid of hard 
spheres with zero dipole moment as a function of the 
packing density 77 = (w/6)p*. For the Carnahan-Starling 
equation of state, on has f{rj) = (477 — 3?? 2 )/(l — i]) 2 . (3 — 
1/T in Eq. Q is given in reduced units and thus the coef- 
ficients A = Na and b depend only on the liquid density 
through two-particle and three-particle perturbation in- 
tegrals, a = (p*/G)I {2) (p*) {L b= (p*/9)j( 3 )(p*)//( 2 ) (/>*), 
tabulated by Larsen et al [25j |. 

Following Freed the density of states can be ob- 
tained by inverse Laplace transformation of Eq. ^ in 
which we use /3F((3) in the Pade form given by Eq. ([4|. 
The inverse Laplace transform is calculated by expand- 
ing exp(— j3F(/3)) in powers of the second summand in 
Eq. (|4|) and performing pole integration. The result can 
be converted to a closed-form equation: 



0(e) = (by/1 + e/c) 1 exp [N(f(r)) - e/b - 2c/b)] 
h (2{c/b)Ny/l + e/c) , 



(5) 



where c = a/b, I\{x) is the first-order modified Bessel 
function, and e = Ea 3 /Nm 2 is the reduced internal 
energy per particle. Since the argument of the Bessel 
function in Eq. ([5]) is proportional to the number of 
particles A, its thermodynamic-limit expansion gives 
the following expression for the enumeration function 




FIG. 1: Enumeration function <r(e) (a) and distribution func- 
tion P(e) (b) of the fluid of dipolar hard spheres at p* = 0.8 
and the temperatures T (from left to right in (b)): 0.125, 
0.167, 0.2, 0.25, 0.33, 0.5, 1.0, 2.0, 10.0. The simulation points 
in (a) are obtained from P(e) within 95% of the maximum 
probability at each temperature. 



<r(e) = A- 1 ln[0(e)]: 



a(e) = f{i 1 )-b- 1 {^T~e-^Y 



(6) 



The enumeration function in Eq. ([6]) is Gaussian near 
its top, e ~ 0, where the dipole-dipole interactions are in- 
significant and the limit of a hard-sphere fluid is reached. 
It deviates from the parabolic shape in its low-energy 
wing, in particular close to the low-energy cutoff at 
eo = — c (Fig. [Iji). Depending on the parameters, Eq. 
© gives two possible resolutions of the entropy crisis 
of low-temperature fluids [l|. When f{rf) — c/b > 0, the 
enumeration function has an infinite derivative at eo and, 
according to Stllinger's arguments @, there is no ideal 
glass transition. When f(rf) — c/b is strictly positive, 
the entropy is non-zero at the cutoff energy, a situation 
observed for network glass formers (26j. Finally, when 
f(r]) — c/b < 0, the ideal glass arrest <j(eic) — happens 
before the cutoff is reached, and the temperature of ideal- 
glass transition Tig, e iG = e (Tic) is positive. The use of 
Carnahan-Starling hard-sphere free energy, and pertur- 
bation integrals from Ref. [2f| gives f(rf) — c/b = —1.19 at 
p* = 0.8 and T IG ~ 0.073. 

The enumeration function from Eq. © is compared 
to the results of MC simulations in Fig. The simu- 
lation points were obtained by patching together central 
parts (within 95% of the maximum) of the distribution 
functions P(e) cx exp(A(er(e) — (3e)) at different temper- 
atures. Since P(e) defines <r(e) up to a constant, this 
comparison only serves to show the non-Gaussian shape 
of the curve. This non-Gaussian form of cr(e) forces the 
width of P(e) to decrease with cooling (Fig. [Tp). This 
behavior is quite distinct from the prediction of the Gaus- 
sian landscape model in which the width is temperature- 
independent [l5l ]. 

The temperature dependence of the width is directly 
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FIG. 2: Average energy e(T) (a), heat capacity cy (b), and 
the configurational entropy s c (c) vs 1/T. The points in (a) 
are the results of MC simulations: open points for the internal 
energies and close points for the energies of inherent struc- 
tures. The solid line in (a) shows the Pade approximation 
of Eq. ([8]). The dashed horizontal lines show the cutoff en- 
ergy eo = —2.07, the energy of non-interacting dipolar chains 
e = —2.4 (middle line), and the energy of the closed-packed 
fee antiferroelectric crystal e^c = —2.56. The inset in (a) 
shows e(T) from simulations, the inset in (b) shows the ne- 
matic order parameter Pi (squares) and the Binder parameter 
Mb (diamonds). The solid line in (c) is the configurational 
entropy s c (T) = cr(e(T)) in the Pade approximation [Eqs. 
([5]) and ©]. The dots show the configurational entropy [Eq. 
©] with e(T) from simulations. The dashed line in (c) is 
the thermodynamic excess entropy over that of the ideal gas 
calculated by temperature integration of simulated cv(T). 



related to the heat capacity since 

(Se) 2 



P(e) oc exp 



2c v T 2 



(7) 



where 8e is the energy fluctuation and cy is the excess 
constant-volume heat capacity over that of the ideal gas. 
The high-temperature portion of the heat capacity, cor- 
responding to the distributions shown in Fig. [Tp, scales 
linearly with the inverse temperature, cy oc 1/T (Fig. 
[2p). This hyperbolic temperature scaling, often docu- 
mented for structural glass formers at constant pressure 
[27| , results in a linear temperature scaling of the squared 
width when the distribution P(e) is fitted to a Gaussian 
function. This behavior was predicted by models of su- 
percooled liquids in terms of configurational excitations 

The calculation of the energies of inherent structures 
by conjugate gradient minimization of configurations 
along simulated trajectories results in energies <j){T) just 
slightly below e(T) (closed points in Fig. Wfr)- This means 



that the long-range dipolar forces produce essentially a 
mean-field potential for the local translations and rota- 
tions and Eq. ^ for the density of internal energies can 
be used for the density of inherent structures as well, 
fl(e) ~ fy(e). 

The density of states from simulations approaches the 
ideal-glass state even steeper than Eq. ^ (Fig.QJi). This 
trend results in a slightly steeper temperature drop of 
e(T) from simulations compared to the result of the Pade 
approximation (Fig. [2^,): 



e(T) = -c+c(l + b/T)- 



(8) 



In particular, instead of leveling off while approaching 
the cutoff energy eo = —2.07, the system undergoes a 
first-order liquid-liquid transition (see below) accompa- 
nied by a discontinuous drop of the internal energy and 
a sharp peak in the heat capacity (Fig. [2k, b). The orien- 
tational structure of the fluid also changes as indicated 
by a step-wise increase in the nematic order parameter 
Pi calculated as the largest eigenvalue of the Q-tensor 



Q = (27V)- 1 ^(3e J e J -I), 



(9) 



where hj is a unit vector along the dipole of particle j. 
The ferroelectric order parameter Pj (normalized total 
dipole moment of the liquid M) remains close to zero in- 
dicating that the liquid remains unpolarized. The Binder 
parameter [H], M B = 1 - (M 4 )/3(M 2 ) 2 , also shows 
downward spikes pointing to first-order phase transitions 
(inset, Fig. [2p). The internal energy at the first spike of 
cy drops from the level eo to the energy approximately 
equal to the energy of non-interacting dipolar chains [29( 
e = —2.40 (middle dashed line in Fig.[3k). This energy is 
slightly above the energy of close-packed antiferroelectric 
fee crystal e^c = —2.56 [3fJ. 

The configurational entropy per particle s c (T) = 
a(e(T)) from Pade thermodynamics [Eqs. © and ©] 
deviates from the hyperbolic functions at low tempera- 
tures approaching the ideal glass transition at a positive 
temperature (Fig. [^t). In contrast, when e(T) from sim- 
ulations is substituted into Eq. ©, the decay of s c (T) is 
faster, reaching the ideal-glass state near the point of the 
first-order transition. The DHS fluids thus loses stabil- 
ity when it runs out of configurations when approaching 
the ideal-glass transition and transforms into a different 
thermodynamic state via a first-order transition. 

The properties of the low-temperature fluid phase are 
peculiar. Both the density (ex (|p/c| 2 )) and polarization 
(oc (|Mfc| 2 )) structure factors do not show crystalline spa- 
tial or orientational order. However, snapshots of simu- 
lated configurations (not shown here) indicate the ex- 
istence of chains of dipoles aligned head-to-tail. These 
chains are persistent for (5 — 20) x I0 3 MC cycles result- 
ing in the overall slow convergence of simulations. The 
chains of parallel dipoles are arranged in bundles with the 
locally body-centered tetragonal (bet) structure in which 
the two chains are displaced relative to each other by the 
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FIG. 3: Pair distribution function of transverse interparti- 
cle separations (relative to the nematic director) at tempera- 
tures indicated in the plot. The high-temperature results for 
T = 0.167 (dashed line) and T = 0.25 almost coincide on the 
scale of the plot. The distribution of distances parallel to the 
director in the inset shows the transition to the smectic phase 
of parallel chains with the local bet order. 

hard-sphere radius. This arrangement has the lowest en- 
ergy for fluids forming columnar phases of parallel dipo- 
lar chains [29|, [3l| . A body-centered orthorombic lattice 
was suggested as the ground state for Stockmayer ferro- 
electrics [23j]. In the present case, all the chain bundles 
are oriented along the same director, but the net moment 
is compensated among the oppositely oriented bundles. 

The columnar structure is well seen in the pair corre- 
lation function of transverse separations, which is very 
sensitive to columnar order [32] (Fig. [3]). At the phase 



transition, a peak at zero transverse projection indicates 
the emergence of the columnar phase. At the same time, 
a smaller peak, present at r/a = 1 in the isotropic dipolar 
fluid, shifts to r/a = v3/2 characteristic of the closest 
bet distance [2{| . The one-dimensional dipolar chains ex- 
perience strong Landau-Peierls fluctuations with the or- 
thogonal mean-square displacement scaling as (rj_) oc IT 
[33[ | . where I is the average length of the chain. The fits 
of the initial portion of g(r±) result is almost invariant 
(rj_) ~ 0.07<7 2 suggesting that the length of the chains 
increases as 1/T. 

The nematic columnar phase transforms into a smec- 
tic phase with further cooling as indicated by the second 
peak of the heat capacity, a downward spike of the Binder 
parameter (Fig. [2t>), and the appearance of clear density 
modulation in the distribution function of molecular sep- 
arations parallel to the director (inset in Fig. [3]). 

In conclusion, a new model of thermodynamics of low- 
temperature fluids offers a description of non-Gaussian 
landscape in a good agreement with simulations. For 
the parameters of the DHS fluids studied here the model 
predicts the excitation profile e(T) to terminate at the 
ideal-glass transition. The DHS fluid nearly avoids this 
point through a first-order transition to a fluid phase with 
columnar order. 

This research was supported by the the Air Force 
(FA9550-06-C-0084) and NSF (CHE-0616646). 



[IT] 



C. A. Angell, K. L. Ngai, G. B. McKenna, and S. W. [18 
Martin, J. Appl. Phys. 88, 3113 (2000). 
I. M. Lifshitz, A. Y. Grosberg, and A. R. Khokhlov, Rev. [19 
Mod. Phys. 50, 683 (1978). 

J. D. Bryngelson and P. G. Wolynes, Proc. Natl. Acad. [20 
Sci. 84, 7524 (1987). 

K. F. Freed, J. Chem. Phys. 119, 5730 (2003). [21 
F. H. Stillinger and T. A. Weber, Phys. Rev. A 25, 978 
(1982). [22 

F. H. Stillinger, J. Chem. Phys. 88, 7818 (1988). 
S. Biichner and A. Heuer, Phys. Rev. E 60, 6507 (1999). [23 
S. Sastry, Nature 409, 164 (2001). 
J. H. Gibbs and E. A. D. Marzio, J. Chem. Phys. 28, 373 [24 
(1958). [25 

G. Adam and J. H. Gibbs, J. Chem. Phys. 43, 139 (1965). 
X. Xia and P. G. Wolynes, Proc. Nat. Acad. Sci. 97, 2990 [26 
(2000). 

A. J. Moreno, I. Saika-Voivod, E. Zaccarelli, E. L. Nave, [27 
S. V. Buldyrev, P. Tartaglia, and F. Sciortino, J. Chem. 
Phys. 124, 204509 (2006). [28 
M. Mezard, G. Parisi, and M. Virasoro, Spin Glass The- 
ory and Beyond (World Scientific, Singapore, 1987). [29 
T. R. Kirkpatrick, D. Thirumalai, and P. G. Wolynes, [30 
Phys. Rev. A 40, 1045 (1989). [31 

B. Derrida, Phys. Rev. B 24, 2613 (1981). 

D. V. Matyushov and C. A. Angell, J. Chem. Phys. 123, [32 
034506 (2005); cond- mat/0612627. [33 
P. G. Debenedetti, F. H. Stillinger, and M. S. Shell, J. 
Phys. Chem. B 107, 14434 (2003). 



C. G. Gray and K. E. Gubbins, Theory of Molecular Liq- 
uids (Clarendon Press, Oxford, 1984). 

J. J. Weis and D. Levesque, Phys. Rev. Lett. 71, 2729 
(1993). 

M. E. van Leeuwen and B. Smit, Phys. Rev. Lett. 71, 
3991 (1993). 

D. Wei and G. N. Patey, Phys. Rev. Lett. 68, 2043 
(1992). 

D. Wei, G. N. Patey, and A. Perera, Phys. Rev. E 47, 
506 (1993). 

G. T. Gao and X. C. Zeng, Phys. Rev. E 61, R2188 
(2000). 

B. Groh and S. Dietrich, Phys. Rev. E 63, 021203 (2001). 
B. Larsen, J. C. Rasaiah, and G. Stell, Mol. Phys. 33, 
987 (1977). 

A. Saksaengwijit, J. Reinisch, and A. Heuer, Phys. Rev. 
Lett. 93, 235701 (2004). 

R. Richert and A. C. Angell, J. Chem. Phys. 108, 9016 
(1998). 

M. S. S. Challa, D. P. Landau, and K. Binder, Phys. Rev. 
B 34, 1841 (1986). 

R. Tao and J. M. Sun, Phys. Rev. Lett. 67, 398 (1991). 
J. M. Luttinger and L. Tisza, Phys. Rev. 70, 954 (1946). 
A. -P. Hynninen and M. Dijkstra, Phys. Rev. Lett. 94, 
138303 (2005). 

D. Wei, Phys. Rev. E 49, 2454 (1994). 

P. G. de Gennes and P. A. Pinkus, Phys. Kondens. Ma- 

terie 11, 189 (1970). 



